function [U,Stress]=Bernoulli3DFEA(Node,Element,Material,force,constrain)
Dof = 6;EDof = Dof*2;
Dofs=Dof*size(Node,1); 
EleCount=size(Element,1); 
[nforce,~]=size(force);[nconstrain,~]=size(constrain);
KKG=sparse(Dofs,Dofs);FFG=sparse(Dofs,1);U=sparse(Dofs,1);
Stress = sparse(EDof,EleCount);
x = Node(:,1);y = Node(:,2);z = Node(:,3);
BarLength=BarsLength(x,y,z,Element);
for iEle =1:EleCount
        n1=Element(iEle,1);n2=Element(iEle,2);
        R=CoordTransform([x(n1) x(n2)],[y(n1) y(n2)],[z(n1) z(n2)],BarLength(iEle));
        ElementStiffnessMatrix = Bernoulli3DElementKe(Material,R,BarLength(iEle),1);
        % 确定节点1和节点2的自由度范围
        n1_dofs = (n1 - 1) * Dof + (1:Dof);
        n2_dofs = (n2 - 1) * Dof + (1:Dof);
        
        % 将节点1和节点2的自由度合并为一个向量
        ElementNodeDOF = [n1_dofs, n2_dofs];
        KKG(ElementNodeDOF,ElementNodeDOF)=KKG(ElementNodeDOF,ElementNodeDOF)+ElementStiffnessMatrix;
end
if ~isempty(force)
    for i=1:nforce
        m=force(i,1);
        n=force(i,2);
        FFG(Dof*(m-1)+n)= force(i,3);
    end
end

for i=1:nconstrain
    n=constrain(i,1); % 节点号
    d=constrain(i,2); % 局部自由度号
    m = (n-1)*Dof + d ; % 全局自由度索引位置
    FFG = FFG - constrain(i,3) * KKG(:,m); % 载荷列阵均减去相应元素
    FFG(m) = constrain(i,3); % 载荷列阵相应位置更改为已知位移
    KKG(:,m)=zeros(Dofs,1); % 刚度矩阵列置零
    KKG(m,:)=zeros(1,Dofs); % 刚度矩阵行置零
    KKG(m,m) = 1.0; % 刚度矩对角元素阵置 1
end


U = KKG\FFG;



end

function BarLength = BarsLength(x,y,z,ele)
   BarLength=zeros(size(ele,1),1);
   for iEle =1: size(ele,1)
      BarLength(iEle,1)=((x(ele(iEle,2))-x(ele(iEle,1)))^2+(y(ele(iEle,2))-...
      y(ele(iEle,1)))^2+(z(ele(iEle,2))-z(ele(iEle,1)))^2)^0.5;
   end
end

function R = CoordTransform(x,y,z,L)
    x1 = x(1);x2 = x(2);y1 = y(1);y2 = y(2);z1 = z(1);z2 = z(2);
    if x1 == x2 && y1 == y2
        if z2 > z1
            Lambda = [0 0 1 ; 0 1 0 ; -1 0 0];
        else
            Lambda = [0 0 -1 ; 0 1 0 ; 1 0 0];
        end
    else
        CXx = (x2-x1)/L;
        CYx = (y2-y1)/L;
        CZx = (z2-z1)/L;
        D = sqrt(CXx*CXx + CYx*CYx);
        CXy = -CYx/D;
        CYy = CXx/D;
        CZy = 0;
        CXz = -CXx*CZx/D;
        CYz = -CYx*CZx/D;
        CZz = D;
        Lambda = [CXx CYx CZx ;CXy CYy CZy ;CXz CYz CZz];
    end
    R = [Lambda zeros(3,9); zeros(3) Lambda zeros(3,6);
    zeros(3,6) Lambda zeros(3);zeros(3,9) Lambda];
end

function Ke = Bernoulli3DElementKe(prop,R,L,icoord)
    E = prop(1);A = prop(2);Iy = prop(3);Iz = prop(4);G = prop(5);J = prop(6); 
    k1 = E*A/L;
    k2 = 12*E*Iz/(L*L*L);
    k3 = 6*E*Iz/(L*L);
    k4 = 4*E*Iz/L;
    k5 = 2*E*Iz/L;
    k6 = 12*E*Iy/(L*L*L);
    k7 = 6*E*Iy/(L*L);
    k8 = 4*E*Iy/L;
    k9 = 2*E*Iy/L;
    k10 = G*J/L;
    ke = [k1 0 0 0 0 0 -k1 0 0 0 0 0;
          0 k2 0 0 0 k3 0 -k2 0 0 0 k3;
          0 0 k6 0 -k7 0 0 0 -k6 0 -k7 0;
          0 0 0 k10 0 0 0 0 0 -k10 0 0;
          0 0 -k7 0 k8 0 0 0 k7 0 k9 0;
          0 k3 0 0 0 k4 0 -k3 0 0 0 k5;
         -k1 0 0 0 0 0 k1 0 0 0 0 0;
          0 -k2 0 0 0 -k3 0 k2 0 0 0 -k3;
          0 0 -k6 0 k7 0 0 0 k6 0 k7 0;
          0 0 0 -k10 0 0 0 0 0 k10 0 0;
          0 0 -k7 0 k9 0 0 0 k7 0 k8 0;
          0 k3 0 0 0 k5 0 -k3 0 0 0 k4];
   
    switch icoord
        case 1
            Ke = R'*ke*R;
        case 2
            Ke = ke;
    end
   
end















